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Abstract 

In this paper, approximate Linear Minimum Variance (LMV) niters for continuous-discrete state space 
models are introduced. Such filters are obtained by means of a recursive approximation to the predictions 
for the first two moments of the state equation. It is shown that the approximate filters converge to the 
exact LMV filter when the error between the predictions and their approximations decreases. As particular 
instance, the order-/3 Local Linearization filters are presented and expounded in detail. Practical algorithms 
are also provided and their performance in simulation is illustrated with various examples. The proposed 
filters are intended for the recurrent practical situation where a nonlinear stochastic system should be 
identified from a reduced number of partial and noisy observations. 

1 Introduction 

The estimation of unobserved states of a continuous stochastic dynamical system from noisy discrete observa- 
tions is of central importance to solve diverse scientific and technological problems. The major contribution 
to the solution of this estimation problem is due to Kalman and Bucy [261 27], who provided a sequential and 
computationally efficient solution to the optimal filtering and prediction problem for linear state space models 
with additive noise. However, the optimal estimation of nonlinear state space models is still a subject of active 
researches. Typically, the solution of optimal filtering problems involves the resolution of evolution equations 
for conditional probabilistic densities, moments or modes, which in general have explicit solutions in few partic- 
ular cases. Therefore, a variety approximations have been developed. Examples of such approximate nonlinear 
filters are the classical ones as the Extended Kalman, the Iterated Extended Kalman, the Gaussian and the 
Modified Gaussian filters [17]; and other relative recent as the Local Linearization [36l[24|, the Projection [3J 
and the Particle filters [TT] methods. 

In a variety of practical situations, the solution of the general optimal filtering problem is dispensable 
since the solution provided by a suboptimal filter is satisfactory. This is the case of the signal filtering and 
detection problems, the system stabilization, and the parameter estimation of nonlinear systems, among others. 
Prominent examples of suboptimal filters are the linear, the quadratic and the polynomial one, which have been 
widely used for the estimation of the state of both, continuous-continuous [301 Ell SO] and discrete-discrete 
[21 [Ml SOI EH U\ models. In the case of continuous-discrete models, exact expressions for Linear Minimum 
Variance filter (LMV) have been derived as well [53], but they are restricted to linear models. For nonlinear 
models, this kind of suboptimal filter has in general no exact solution since the first two conditional moments of 
the state equation has no explicit solution. Therefore, adequate approximations are required in this situation 
as well. 
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In this paper, approximate LMV filters for continuous-discrete state space models are introduced and their 
convergence is studied. The filters are obtained by means of a recursive approximation to the predictions 
for the first two moments of the state equation. Based on the well-known Local Linear approximations for 
the state equation, the order-/3 Local Linearization filters are presented as a particular instance. Their con- 
vergence, practical algorithms and performance in simulations are also considered. It is demonstrated that 
these approximate filters enhance the identification of nonlinear stochastic systems given a reduced number of 
partial and noisy observations, which is a typical situation in practical control engineering. Moreover, they are 
crucial for improving the parameter estimation of this type of state space models by means of the innovation 
method [37J 123 HI] which, defined in term of the exact LMV filters, is a prediction estimation method [53] 
with convenient statistical properties. 

The paper is organized as follows. In section 2, basic notations and results on LMV filters, Local Linear 
approximations and Local Linearization filters are presented. The general class of approximate LMV filters 
is introduced in section 3 and its convergence stated. In section 4, the order-/? Local Linearization filters are 
presented and their convergence analyzed. In the last two sections, practical algorithms for these filters and 
their performance in simulations are considered. 

2 Notations and Preliminaries 

Let (JljJ 7 , P) be a complete probability space, and {J^t, t > to} be an increasing right continuous family of 
complete sub cr-algebras of T . Consider the state space model defined by the continuous state equation 

ni 

<fjt(t) = f(t,x(i))di + ^g i (t,x(i))dw i (i), forte[i ,T], (1) 

i=l 

and the discrete observation equation 

z tk = Cx(i fe )+e tfe , forfc = 0,l,..,iV, (2) 

where f, g^ : [to,? 1 ] x R d — > M. d are functions, w = (w 1 , . . . , w m ) is an m-dimensional J^-adapted standard 
Wiener process, {e tfc : e tk ~ Af(0, X), k — Q,..,N} is a sequence of r-dimensional i.i.d. Gaussian random 
vectors independent of w, and C an r x d matrix. Here, it is assumed that the N time instants tk define an 
increased sequence {£}jv = {tk ■ tk < tk+x, tpj = T, k = 0,1, .., N — 1}. Conditions for the existence and 
uniqueness of a strong solution of ([I]) with bounded moments are assumed. 

Let x t / p — E(x(t)/Z p ) and Qt/ P = E '(x(i)x T (t) / Z p ) be the first two conditional moments of x with p <t, 
where E{.) denotes the mathematical expectation value, and Z p = {z, tk : t k < p} is a time series with 
observations from ph. Further, let us denote by 

u t/p = £((x(t) - x t/p )( x (i) - * t/p y/z p ) 

= Qt/p - x t/p x J/ p 

the conditional variance of x. 

Denote by Cp(M d ,K) the space of I time continuously differentiable functions g : M. d — > K for which g and 
all its partial derivatives up to order I have polynomial growth. 

2.1 Linear minimum variance filtering problem 

According to [2] |44] HH [IT] the linear minimum variance filter x t k+1 /t fc+I for a state space model with discrete 
observation ^ is defined as 

x t k+1 /t k+1 = x t fc+1 /t fc + G t(c+1 (z (fc+i - Cx tk+1 / tk ), 
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where the filter gain G tk+1 is to be determined so as to minimize the error variance 

£((x(t fc+ i) - x tfc+l/tjfc+1 )(x(t fe+1 ) - x tfc+l/tfc+1 ) T ). 
This yields to the following definition. 

Definition 1 The Linear Minimum Variance filter for the state space model is defined, between obser- 

vations, by 

' /X " ' - E(f(t,x)/Z t ) (3) 



dt 

dJJ m 

= £(xF(t, x)/Z t ) - x t/t £(F(t, x)/%) + E(f (t, x)xT/Z t ) - £(f (t, x)/Z t )xJ /t - x)gT(i, x)/Z t ) 



(4) 

/or i £ (tk,tk+i), and by 

x t k+1 /t k+1 = x t fc+ i/t fc + G tfc+1 (z tk+i - Cx (t+1 / t J (5) 

u u-+i/t fc +i = u t fc +i/t fc - G tfc+1 CU tfc+l/tfc (6) 
/or eac/i observation at tk+i, with filter gain 

Gt k+1 = U tfe+l/tfc CT(CU tfc+l/tfc CT + S^J" 1 (7) 

for all tk,tk+i € {£}jv- The predictions x t / tk , Ut/t k are accomplished, respectively, via expressions 
initial conditions x tfc / tfe and U t( ./ t( . /or £ € (ifc,ifc+i] and £fc ■ ifc+i G Ww- 

2.2 Local Linearization filters 

A key component for constructing the Local Linearization (LL) filters is the concept of Weak Local Linear 
(WLL) approximation for Stochastic Differential Equations (SDEs) [Ti?I |2~1] . 

Let us consider the SDE (Jfj), and the time discretization {r) h — {r„ : n = 0, 1, . . . , M} of the time interval 
[to , T] with maximum stepsize h defined as a sequence of times that satisfy the conditions to = tq < t\ < ■ ■ ■ < 
tm = T, and max(r n+ i — r„) < h < 1 for n = 0, . . . , M — 1. Further, let 

n 

n t = max{n = 0, 1, . . . , M : r„ < t and r„ £ (t)^} 

for all i G [t , T]. 

Definition 2 for a given time discretization (r) h , the stochastic process y = {y(t), t € [to,T]} is called order- 
(3 (= 1,2) Weak Local Linear approximation of the solution of |IJ) if it is the weak solution of the piecewise 
linear equation 

m 

dy(t) = (A(T„ t )y(t)+a^(t;T„ t ))dt + ^(B i (T„ t )y(t)+bf(t;r nt ))dw i (t) (8) 

i=l 

for all t £ (t„,t„ + i] and initial value y(io) = x (io); where the matrices A,Bj are defined as 

dt(s,y(s)) dgj(s,y(s)) 
A[s) = — and B i (s) = — , 



and the vectors sP , b 



' i 



as 



f (*, y(s)) - ^%(«) + »?^fJi(t _ a) for = 1 

^(t;s) + l t [G( S ,y( S ))GT( S) y( s ))p Sp ^0 1 {t-s) /or/3 = 2 
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and 

f ft(»,y(«)) - » y(*) + 9eds d * {s)) (t- s ) for = 1 

b ^ (<;s)= bHt;«) + | t [G( S; y( s ))GT( s , y ( s ))p g!|M (t _ s) /or/3 = 2 

/or a// s < t. Here, G — [g x , . . . , g m ] is an d x m matrix function. 

The drift and diffusion coefficients of the equation {[sj) are, respectively, weak approximations of order j3 to 
the drift and diffusion coefficients of the equation (JlJ obtained from the Ito- Taylor expansion of order (3. That 



sup 

s<t<s+h 



\E(g(f(t, y(t))) - E(g(A(s)y(t) + a"(*; s)))\ < Ch? 



and 



sup 

s<t<s+h 



E(g(gi(t, y(t))) - E(g(B t (s)y(t) + bf (t; a))) 



<Ch? 



for all h > and s G [toj^" — h], where g £ Cp' 3+1 ' ) (]R d ,M) and C is a positive constant. 

Explicit formulas for the conditional mean y t i p and variance V t i p of y were initially given in |231 124) and 
simplified later in fT8] , 

The conventional Local Linearization filters for the model ([l])-([2]) are obtained in two steps [23]: 1) by 
approximating the nonlinear state equation by the Local Linear approximation ^ on a time partition (r) h = 
{t}N] and 2) by applying the linear minimum variance filter [23] to the resulting linear continuous-discrete 
model. This yields to the following. 

Definition 3 Given a time discretization (r) h = {£}at, the Local Linearization filter for the state space model 
|Ip-Q) is defined, between observations, by 

d ^=A(t nt )y t/t + ^(t;t nt ) (9) 
= A(t nt )V t/t + V t/t AT(t nt ) + ^Bi(i nt )V t/t Bj(t„ t )+B(t;t nt ) (10) 

i— 1 

for all t E (tk,t k+1 ), and by 

Yt k+1 /t k+1 =Yt k+1 /t k +K tfc+1 (z tfc+i -Cy tk+l/tk ) (11) 

Y t k+1 /t k+1 = Y t k+1 /t k ~ K tk+1 CV tk+l/tk (12) 
for each observation at tk+i, with filter gain 

K tk+1 = V tfc+l/tfc CT(CV tfe+l/tfc CT + Etfc+i )-i (13) 



for alltk,tk+\ S {t}N- The predictions y t / tk andV t / tk are accomplished, respectively, via expressions ^-{10) 
with initial conditions yt k /t k an d ^t k /t h f or t € (*fc> tk+i\- Here, 



B(t; s)=fl B,( S )y t/t yJ /t Bj( S ) + B % ( S )y t/t (b^t; a ))T + bf (t; s)yJ /t B](t) + b,(i; S )(bf (t; s)) T , 

i=l 

and £/ie matrices A, B^ and </ie vectors a, bf are defined as in the WLL approximation |#|). 
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Both, the Local Linear approximations and the Local Linearization filters have had a number of important 
applications. The first ones, in addition to the filtering problems, have been used for the derivation of effective 
integration [THl G3 [51 [52] and inference [55J 251 QUI EH [TFJ methods for SDEs, in the estimation of distribution 
functions in Monte Carlo Markov Chain methods [391 [331 HI and the simulation of likelihood functions [3SJ . 
The second ones have played a crucial role in the practical implementation of innovation estimators for the 
identification of continuous-discrete state space models |38[ 125] . In a variety of applications, the innovation 
methods have shown high effectiveness and efficiency for the estimation of unobserved components and unknown 
parameters of SDEs given a set of discrete observations. Remarkable is the identification, from actual data, of 
neurophysiological, financial and molecular models among others (see, e.g., }4"1 [51 |2"TI liTl WI\ ) . 

3 Approximate Linear Minimum Variance filters 

Let (t)l be a time discretization of [id? 1 ] sucn that (r) h D {t}N, and y n the approximate value of x(i n ) 
obtained from a discretization of the equation |TJ) for all r n 6 (j) h . Let us consider the continuous time 
approximation y = {y(t), t £ [to,T] : y(T„) = y„ for all t„ 6 (t)/,} of x with initial conditions 



E (y(t )\F t0 ) = E \x(to)\ Ttoj and E (y(t )y T (i ) |^t J = E (^x(£ )x T (t ) | T ta 
satisfying the bound condition 

E^y(t)\ 2<! \^ \ <M X (14) 
for all t € [toi? 1 ]; and the weak convergence criteria 



sup 



E[g(x(t)) \F tk - E \g{y{t)) LF tj 



< M 2 ^ (15) 



for all fc = 0, 1, .., N — 1, where g € Cp /3+1 ' (R d , R), Mi, M 2 are positive constants, and q = 1, 2.... The process 
y defined in this way is typically called order-/? approximation to x in weak sense [28] . 

When an order-/? approximation is used for the solution of the state equation 0, the following approximate 
filters can be naturally defined. 

Definition 4 Given a time discretization (j) h D {i}jv ; the order-fi Linear Minimum Variance filter for the 
state space model |I])-|IJ) is defined, between observations, by 

y t/t = E(y(t)/Z t ) and V t/t = £(y(t)y T (t)/Z t ) - y t/t yj /t (16) 
for all t £ (ife, tk+i), and by 

Yt k+1 /t k+1 =Yt k+1 /t k + K t fc+ i( z t fc+1 - Cy u+l/t J, (17) 

V tfe+1 / tfc+1 = V tk+1 /t k - K tfc+1 CV tfc+l/tfc , (18) 
for each observation at tk+i, with filter gain 

K tk+1 = V tk+l/tk C\CV tk+l/tk Ci + S^J" 1 (19) 

for all t k ,t k+ i G {t} N . The predictions y t/tk = E(y(t)/Z tk ) and V t/tk = E{y{t)yt (t) / Z tk ) - y f/tfc y t T /tfc , ic&SA 
initial conditions y tk /t h an d V tfc / tfc7 are defined for all t 6 (ifc,ifc + i] andtk,tk+\ G {<}v. 

The next result deals with the convergence rate of the above approximate filters to the LMV filter. 
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Theorem 5 Letx t / p andU t / p be the exact LMV estimates defined in Then, between observations, the 

approximate order-/3 LMV estimates satisfy 



|x t/t - y t / t \ < Kh fj and |U t/t - V t/t \ < KhP 
for all t G (tk,tk+i) and, at each observation t^+i, 

| X *fc-t-i/tfc+i y^fc + i/ifc + i 

for all tfc,tfc+i G {£}jv 7 where K is a positive constant. For the predictions 



< Eh? and |U tfc+l/tfc+1 - V tk+l/tk+1 \ < Kh? 



|x t/4fc - y t/tfc | < KhP and |U t/tfc - V t/tk \ < Kh? 



(20) 



(21) 



(22) 



hold for all t G (tfc, tfc+i] STid t^, tk+x G {i}./v- 



Proof. Let us start proving inequalities (22| and (20). Condition (151 directly implies |x f / tfc — Yt/t k \ < A^/i' 3 
and | Qt/*„ - Pt/t fe | < M 2 ft^, where P t/tk = £(y(i)y T (i)/Z t J. Condition {li} implies that |y t / tfe | and |P t /tJ 
are bound. From the definition of variance in terms of the first two moments follows that 



|U t/tfc - V t/tk I < \Q t /t h ~ p t/t k I 



x «/^ x I/* fc _ yt/t k y]/ tk 



Since 



x t/t fe x J/ t|! - y t /t k yj/ tk 



< 



x */tfc x t/* fc ~ *t/t k y}/ tk +^t/t k yj/ tk - yt/t k y] /tk 
x */tf=( x t/t fc ~ y J t/t k ) + ( x */tf= _ yt/t k )y] /u 



< (l x */tJ + |y*/* fe ] ) l x t/* fc - y*/*kl 



and \xt/t k I is a ls° bound, there exists a positive constant L such that |U t / tfc — V t / tfc | < Lft,' 9 . Inequalities (20) 
can be derived in the same way. 

For the remainder inequalities follows this. From ^ and (17), it is obtained 



l x *fc+i/tfc+i ~ ytfc+i/tfc+i I - l x *fc+i/t fc - yt k+1 /t k I + G ^+i( z t fc+1 _ Cx t k+1 /t k ) - K t fc+1 ( z t fc+1 - c yt k+1 /t k ) 

< (! + | G t fc+1 c |) | x t fc+1 /t fc ~ yt k+1 /t k I + (|z ffc+1 1 + Cy tk+i/tk ) \G tk+1 - K tfc+1 1 . 
From ^ and (18), it is obtained 

\^t k+1 /t k+1 - V tfe+1( / ijb +1 J < |U tfc+1 / tfc - V tfc+1 / tfe | + |G tfc+1 CU tfc+1 / tfc - K tfc+1 CV tfc+1 / tfc | 

< (1 + | G t fc+ i C |) |U tfc+1 / tfc - V tfc+l/tfc | + |CV tfe+1 / t J |G th+1 - K tfc+1 1 . 

By rewriting ( 19 ) and (JtJ) as 

K tfe+1 (CV tfc+l/tfc CT + S tfe+1 ) = V tfc+l/tfc CT 

and 

G tfc+1 (CU tfc+l/tfe CT + E tk+1 ) - G tfc+1 (CV tfc+l/tfc CT + S u+1 ) + G ifc+1 (CV ffc+l/tfc CT + S tfc+1 ) = U tfc+l/tfc CT, 
and subtracting the first expression to the second one follows that 

S t fc +J - G t k+ i C (~V tk+l/tk - U t(c+l/t JC T + (V tk+1 / th - Vt k+1 /t k 



Thus, 

(G tk+1 -K tk+1 ) = (I-G th+1 C)(V th+i/th - V tfc+l/tfc )C T (CV tt+l/tfc CT +S tfc+1 )- 1 

and 

\G tk+1 -K tk+1 \ < |(I-G tfc+1 C)| |CT(CV tfe+l/tfc CT +53 tfc+1 )- 1 | \V tk+l/tk - V th+l/tk \ . 

From the above inequalities, and taking into account that | V t / tfc | , Gt fc+1 , |S| and |C| are also bound, 
there exists a positive constant M such that |x ffc+l/tfc+1 - y tfc+l/tfe+1 | < Mh^ and |U tfc+l/tfe+1 - V tk+l/tk+1 \ < 
MhP . This concludes the prove. ■ 

The result of the Theorem above has been stated for models with linear observation equation. However, by 
following the procedure proposed in [35], it can be easily applied as well to models with nonlinear observation 
equation. 

For illustrate this, let us consider the state space model defined by the continuous state equation ([I]) and 
the discrete observation equation 

z tk = h(t k , x(t fc )) + e tk , for k = 0, 1, .., N, (23) 
where e tk is defined as in ^ and h : K x R d — > W is a twice differentiable function. By using the Ito formula, 

fc=l s=l k,l=l s=l 1=1 

m 

= p> dt + ^aidv/ s 

8=1 

with j = l,..,r. Hence, the state space model ([T])+(|2"3"|) is transformed to the following higher-dimensional 
state space model with linear observation 

m 

dv(t) = a(t, v(t))dt + ^2 h dt, v(t))dw l (t), 

i=l 

z tk =Cv(tjfc) + e tfc , for k = 0,1,.., N, 

where 





X 




' f " 


, b, = 


Si 


V = 












h 




. P _ 





and the matrix C is such that h(i/-, x(tfe)) = Cv(tfc). 

In this way, the state space model (flj) + (p3| is transformed to the form of the state space model Q-Q, 
and so Theorem [5] can be applied. 

4 Order-/3 Local Linearization niters 

In principle, according to Theorem [5j any kind of approximation y converging to x in a weak sense can be 
used to construct approximate LMV filters (e.g., those in [2E]). Therefore, additional selection criterions could 
be taking into account for this purpose. For instance, high order of convergence, efficient algorithm for the 
computation of the moments, and so on. In this paper, we elected the Local Linear approximation |8]) for the 
following reasons: 1) its first two conditional moments have simple explicit formulas that can be computed by 
means of efficient algorithm (including high dimensional state equations) [23 , [231 HE] ; 2) its first two conditional 
moments are exact for linear state equations in all the possible variants (with additive and/or multiplicative 
noise, autonomous or not) |23j : and 3) it has an adequate order of weak convergence for state equations with 
additive noise [6]. 
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Once the order-/? Local Linear approximation Q is chosen for approximating the state equation (JlJ, the 
well know ordinary differential equations for the first two moments of linear SDEs [T] can be directly used to 
define the following filter. 

Definition 6 Given a time discretization (r) h D {t}N, the order-fi Local Linearization filter for the state space 
model is defined, between observations, by 



dy t /t 
dt 



= A(T„ t )y t/t +a' 9 (i;T„ t ) 



dP m 

-£i = A(r Bt )P t/t + P i/t AT(r„ t )+^B i (r n JP t/t BT(r„J + B(t;r nt ) 



v t/i = P t /t ~ yt/tyl/t 

for all t e {t k .t k+ i), and by 

yt fc+ i/t fc +i =yt k+ i/t k + K t fc+ i( z tfc+1 - c yt k+1 /t k ) 

Vt k+ i/t h+1 = Vt k+1 /t k ~ K t)c+1 CV tfc+1 / tfc 
for each observation at tk+i, with filter gain 

K tk+1 = V tfc+l/tfc CT(CV tfc+l/tfc CT + E^)" 1 



(24) 

(25) 
(26) 

(27) 
(28) 

(29) 



for all tk,tk+i € \t\ 7v • The predictions y t / tk , ~Pt/t k ® n d'Vt/t k accomplished, respectively, via expressions 



24)-(26) with initial conditions yt k /t k an d Pt k /t k for t € (ifc,ifc+i] an d tk,tk+i € {£}jv- Here, 



B(t; s) = a"(t; S )yJ /t + y t/t (a^t; s)) T + £ B,( S )y 4/t (bf (t; fl ))T + bf (t; S )y t T /t Bj( S ) + bf (t; S )(bf (t; S )) T , (30) 

i=i 

and the matrix functions A,B; and the vector functions a^,b^ are defined as in the WLL approximation |$|). 

The approximate LL filter ( 24 1-( 29 ) reduces to the conventional LL filter ([9])-( 13 1 when [r) h = {t}^. For 
linear state equations with multiplicative noise, the LL filter (24)- (29 1 reduces to LMV filter proposed in [23] . 
whereas for linear state equations with additive noise, the LL filter (24)- (29 1 reduces to the classical Kalman 
filter. 

According with Theorem [5] the approximate LL filter (24|-(29| will inherit the order of convergence of 
the WLL approximation As it was mention before, the weak convergence rate of that approximation was 
early stated in 6j for SDEs with additive noise. For equations with multiplicative noise, this subject will be 
considered in what follows. 

Lemma 7 Suppose that the drift and diffusion coefficients of the SDE |Ip satisfy the following conditions 

f fe ,g£ G C 2 p (l3+1) ([t Q ,T} x R d ,R) (31) 



|f( S) u)|+£(| gi ( S ,u)|+ |g?(s,u)g<( S ,u)|^)<X(l + |u| 



dt(s,u) 



dt 



k,i=i 



df (s, u) 



<9> 



<9 2 f(s,u) 



9x 2 



6}<K 



(32) 
(33) 
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and 



dgi(s,u) 




<9g 4 (s,u) 




d 2 gi(s,u) 


dt 


+ 


<9x 


+ 


chc 2 



6*<K 



(34) 



for all sG [fo,T], u e R d , andi = 1, ..,m, where K is a positive constant. Then the order- /3 WLL approximation 
|^) satisfies 

(35) 



£ sup |y(t)r Ji <C(l + |y n 

\t <t<T J 

for each q = 1, 2, ... , where C is positive constant. 

Proof. Let us denote the drift and diffusion coefficients of the SDE Q by 

p(i, y (t) ;r nt ) = A(r„ t )y (t) + (i; r„ t ) 



and 



qi(*,y(*);T-n t ) = B,(T„ t ) y (t) + bf (t;r„ t ), 



respectively. 

For each g, the Ito formula applied to |y(i)| 2? implies that 



z t 

|y(t)| 29 = |y(r, it )| 29 + / 2 q \y( S )\ 2 "- 2 y^s)p( S Ms);r nt )d S + Y, / 2<Z |y(*)| 2 ^ 2 y T ( S )q,:( S , yOO^Jdw^s) 

t t 

+ E J q\y(s)\ 2q - 2 \<ii(s,y(sy,r nt )\ 2 ds + Y^ J M<1 1) \y(s)\ 2q - A |y T ( S ) qi ( S ,y( S );r nt )| 2 ds 



for all i € [r„ f ,r„ t+ i]. 

By recursive application of the expression above it is obtained that 



|y(i)| 2 " = |y(*o)| 29 + / 2< Z |y( S )| 29 - 2 yT( s )p( S! y( s );r„Jd S + E / 2g |y( S )| 2? - 2 y"(s) qi (s, y(a);r n .)dw'( a ) 

£ i 
m « m. r. 



to 



- Co 



for all te[t ,T]. 

Theorem 4.5.4 in 28] implies that E [\y{t)\^ q ^j < oo for to < t < T. Hence, the function r defined as 

r(t) = for < t < t and as r(t) = |y(t)| 29 ~ 2 y T (i)q 4 (t, y(£);T„ t ) for < < t < T belongs to the class L\ of 
function C x T— measurable. Then, Lemma 3.2.2 in [28] implies that 



E | J |y(s)| 2 ^ 2 y T ( S )q 4 ( S ,y( S );^JdwX s ) | = o 

\t 
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for all i = 1, ..,m. From this and the previous expression for \y{t)\ 2q follows that 

t 

e( sup |y( M )| 2 «|j- in ) <|y(i )| 2 « + 2g / £ f |y(.s)| 2 «- 2 |yT( S )p( S , y( S ); r„J| | _F to ) ds 

\t <u<t J J V / 

to 

+ «E /^(|y(s)| 29 - 2 |qi( S) y( S );r„J| 2 |j-Jd S 
t 

+ 2g(g-l)f; |s^|y( S )| 29 - 4 |yT( s )q 4 ( s ,y( s );r„J| 2 |j- t0 ^ S . 

i=1 *o 



From conditions ( 32 1-( 34 ) follows that 

|p(*,y(*);r„.)| < K{\y{s)\ + |y(r B .)|) + K P {\ + \y(s)\) +K 

and 



| qi (s,y(s);r n J| < K(\y(s)\ + |y(r„.)|) + i^(l + |y(s)|) + if, 

where 



K for /3 = 1 

if (1 + i-fT) for = 2 

Thus, there exists a positive constant C such that 



|yT( s )p( s ,y( s );r„J| < C(l + \y(s)\ 2 ) + C(l + |y(r„J|^), 

|q,(5,y( S );r„J| 2 < C(l + |y( s )| 2 ) + C(l + |y(r„ s )| 2 ), 
|y T ( S ) qi ( S ,y( S );r„J| 2 < C|y( S )| 2 (l + |y(,s)| 2 ) + C|y( S )| 2 (l + |y(r n J| 2 ), 



and so 



E ( sup \y(u)\ 2q \T t( ) < \y \ 2q + L f E ( sup (1 + |y(«)| 2 ) \y(u)\ 2q - 2 \j"t ) ds, 

\t <u<t / J \t <u<s / 

to 

where L = 2qC(2 + 2qm - m). From the inequality (1 + a 2 )a 2q ~ 2 < 1 + 2a 2<z , 

t 

e( sup |y( U )| 2,? |j" i0 ) < |y | 2 « + L(i-* ) + 2L f £ ( sup |y(u)| 2 « | F to ) ds. 

\t <u<t J J \t <u<s J 

to 

From this and the Gronwall Lemma, the assertion of the Theorem is obtained. ■ 

In what follows, additional notations and results of [28j will be used. Briefly recall us that M. denotes the set 
of all the multi- indexes a = (jx, . . . , ju a )) with j$ € {0, 1, ... , m} and i = 1, . . . , 1(a), where m is the dimension 
of w in ([!]). 1(a) denotes the length of the multi-index a and n(a) the number of its zero components, —a and 
a— are the multi-indexes in M. obtained by deleting the first and the last component of a, respectively. The 
multi-index of length zero will be denoted by v. Further, 

F) d r) 1 d m fi 2 



k=l k,l=lj=l 
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denotes the diffusion operator for the SDE (JlJ, and 

d q 
k = l 

for j = 1, . . . , m. 

Let us consider the hierarchical set 

Tfi = {aeM: 1(a) < {3} 
with f3 = 1, 2; and ^(r^) — {a G M\Tp : -a G Tp} the remainder set of Tp. 

Lemma 8 Let y be the order-fj WLL approximation |5|) ; and z = {z(t), t € [io,? 1 ]} fee t/ie stochastic process 
defined by 

Z(t)=y nt + ^2 I a[^a(T nt ,y nt ;T nt )] Tn ^t+ ^ I a[ A a(^y-'-, T n t )}r nt ,t, (36) 
aeTp/{v} aeB{Tp) 

where, for any given (r„ t ,y„ t ), 

, , J ^...i»'«w-V(a.v;r Bt ) ifji( a) =0 

QlS ' V ' TnJ | L^...L^)-qf i(Q) ( S ,v;r nt ) »/ j I(a) ^ 

is a function of s and v, rnf/i 

p /3 (s,v;r„ t ) = A(r„ f )v + a /3 (s;r„ t ) and qf(s,v;r ra J = B<(r nt )v + bf (s;r„ t ), 

/or a// s G [to)^ 1 ] anc ^ v G ^ > an d matrix functions A, and vector functions a^,bf defined as in the WLL 
approximation Q). Then 

E(g(y(t)))=E(g(z(t))), (37) 
S( 5 (y(t)-y(r nt ))) = S( 5 (z(t)-z(r nt ))) (38) 
for all t € [i ,T] and 5 G Cp (/3+1) (K d , M); and 

la 

/or a^/ a G T^/jV} and t G [io,T], where A Q denotes the Ito coefficient function corresponding to the SDE |Ip. 

Proof. The identities (37)-([38| trivially hold, since (36) is the order- j3 weak Ito- Taylor expansion of the 
solution of the piecewise linear equation Q with initial value y(to) — x (^o)- 

By simple calculations it is obtained that Ito coefficient functions A a corresponding to the SDE ([!]) arc 

\ k e k 

A (0) - 1 ) 



\k „k 



A (i) - ^ • 



A 



a k d c, k , d m c,2 k 



(o,i) dt ^ <9x* 2 ^ ^ dx^x 1 ' 



A fc -YV^ 

2=1 

df k d df k 1 d m <9 2 f fc 

i=l i,i=lj=l 

1=1 
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for a £ F 2 . By taking into account that p' 3 (s,v;r„) and qf(s,v;r„) are linear functions of s and v, it is not 
difficult to obtain that A Q (t„ s , y„ s ; r„ s ) = A Q (r„ s , y ns )t n ,s for all a £ r 2 , which implies (39). ■ 

Note that, the stochastic process z defined in the previous lemma is solution of the piecewise linear SDE Q 
and A Q denotes the fto coefficient functions corresponding to that equation. Therefore, (36) is the Ito-Taylor 
expansion of the Local Linear approximation ([8]). 

The main convergence result for the WLL approximations is them stated in the following theorem. 

Theorem 9 Let x be the solution of the SDE and y the order-fj weak Local Linear approximation of x 
defined by Suppose that the drift and diffusion coefficients of the SDE satisfy the conditions (31)-{34). 
Further, suppose that the initial values of x satisfies 



£(|x H <co 



(40) 



for q = 1,2,..., some constant Co > and all g £ Cp /3+1 ' l (IR d , M). Then there exits a positive constant C g 
such that 

\E(g(x(T)))-E(g(y(T)))\<C g (T-t )hP. (41) 
Proof. For I = 1, 2, . . . , let P, = {p £ {1, . . . , d} 1 }, and let F p : R d ->• K be the function defined as 

/ 

F p (x)=nx p % 



where p = (pi, . . . ,p{] £ Pi. 

By applying Lemma 5.11.7 in [25J to ([8| and taking into account that (36) is the order-/3 weak Ito-Taylor 
expansion of the solution of (fsl) it is obtained 



E ( F p (y„ +1 - y„) - F p ( ^ I a[^ a (T n ,y n ;T n 



<K(l + \y n \ 2r )(T n+1 -T n )h?, 



for all p £ Pi and I = 1, . . . , 2/3 + 1, some K > and r <E {1,2,.. .}, where A Q denotes the Ito coefficient 
function corresponding to ([8]), and h„ = T n+ \ — r n . Further, Lemma [8] implies that 



E F p ( Ia[K{T n ,y n )] TntTn+1 )\F Tn ) \=E Fp( J a[-M T ™> y« - T «)k,r„ +1 ) | J"r n 

\ aei>/M / \ aer>/M 

where A Q denotes the Ito coefficient function corresponding to ([!]). Hence, 



( F p (y„ +1 - y„) - F p ( ^ / a[AQ(r rl ,y n )] TiiiT7i+1 ) | J 7 ^ 



<K(l + \y n \ 2r )(r n+1 -r n )h^ n 



< K(l + max \y k \ 2r )(T n+1 - T n )h p n . 

0<k<n 



On the other hand, Theorem 4.5.4 in [28] applied to Q and Lemma [7] imply 



E ( |y„+i - y n \ 2q \T Tn ) < L{\ + max |y fc | a *)(7- n+ i - r„) 9 

0<A;<7i 



1 2 ^ 



for all < n < M - 1, and 



E[ max |y fe | 29 |J- t0 ) <C(l + |yor"), 
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respectively, where C and L are positive constants. The proof concludes by using Theorem 14.5.2 in [28j with 
the last three inequalities. ■ 

For state equations with additive noise, the order of weak convergence of the WLL approximations provided 
by this Theorem matches with that early obtained in [5J. 

Theorem [9] provides the global order of weak convergence for the WLL approximations at the time t = T . 



Notice further that inequality (41 1 implies that the uniform bound 



sup \E(g(x(t)))-E(g(y(t)))\<C g he 
te[t ,T] 

also holds for the order-/? WLL approximation y, since the order condition provided by Theorem [9] also holds 
uniformly in [to,T] (see Theorem 14.5.1 and Exercise 14.5.3 in [28] for details). 

Convergence in Theorem [9] has been proved under the assumption of continuity for f and gj . If that is 
not the case, the consistency of the WLL discretization has been proved in [5U]. In other practical situations, 
it is important to integrate SDEs with nonglobal Lipschitz coefficients on M. d [3T]. Typically, for such type 
of equations, the conventional weak integrators display explosive values for some realizations. In such a case, 
if each numerical realization of an order-/? scheme leaving a sufficient large sphere 1Z C K d is rejected, then 
Theorem 2.3 in [31] ensures that the accuracy of the scheme is e + O(h^), where e can be made arbitrary small 
with increasing the sphere radius. Therefore, such procedure can be applied to WLL approximation as well. 

Finally, the rate of convergence of the approximate Local Linearization filter is states as follows. 



Theorem 10 Under the assumptions of Theorem^ the approximate order-/3 LL filter (24-)-(29) converge to 
the exact LMV filter (^-{^ with order /?. 

Proof. It is a straightforward consequence of Theorems [9] and [5j ■ 



5 Practical Algorithms 



This section deals with practical implementation of the order-/? LL filter ( 24 )-( 29 ) . Explicit formulas for the 
predictions yt/t k an d P t / tk , an adaptive strategy for the construction of an adequate time discretization {r) h , 
and the resulting adaptive LL filter algorithm are given. 



5.1 Formulas for the predictions 

Let us define the vectors a (Y), ai(r), bi,o(i") and b i) i(r) satisfying the expressions 

aP(t]T nt ) = a (r„ t ) + ai(r„ t )(f - r„ t ) and bf (t;r„ t ) = b ij0 (T„ t ) + h iA (T nt )(t - r„ t ) 
for a given /?. By simplicity, the supraindex j3 is omitted. 



According Theorem 3.1 in [18] . the solution of the differential equations ( 24 1-( 25 ) for the predictions can 
be computed as 



nt — 1 



yt/t k 



y tk/tk + ]T L 2 e M ^)^— ")u T „ A +L 2 e M ^)( t -^)u r „ t , 



(42) 



and 



rit — 1 



vec(P t/th )= ]T L ie M ^)^ + -^)u T „, ffc +L ie M (^)( t --)u r „ t , tfe 

n=n t k 



(43) 
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for all t € [tk,tk+i] and tk,tk+i S {OjVj where the vector u Tttk and the matrices M(r), Li, L2 are defined as 



M(r) = 



^(t) S 5 (r) 64 (r) B 3 (r) S 2 (r) Si(r) 

C(r) I rf+2 

C(t) 










2 









1 





and 



Ll = [ Id 2 °rf 2 x(2d+7) ] > 

in terms of the matrices and vectors 





uec(P r/t J 









r 


1 U T,t k — 












1 


®dx(d 2 +d+2) 


Irf 0dx5 ] 



p (d 2 +2d+7) 



A(t) = A(t)©A(t) + ^ B i (r)0BT( r ), 



C(r) 



A(r) ai(r) A(r)y T , tfc + a (r) 
1 




»(d+2)x(d+2) 



r T = [ lx(d+1) 1 ] 



Bx(t) = vecfatr)) + fa(T)y r/th , B 2 {t) = vec(j3 2 (r)) + fa(r)y T , th , B 3 (r) = «ec(/3 3 (r)), B 4 (r) = /3 4 (r)L and 
B 5 (r)=^(r)Lwith 



Mr 



m 

X^KoWbT^r) +b M (r)bT ( T ) 

m 

^KxWbJ^r) 

m 

a (r) © a (r) + ^ b ii0 (r) (gi B;(t) + B<(r) ® b;, (r) 

m 

ai(r) © ai(r) + J^bi.ifr) ® B;(r) + B,(t) ® b M (r) 



and L = [ Id 0^x2 ] . The symbols vec, © and <S> denote the vectorization operator, the Kronecker sum and 
product, respectively. 

Alternatively, see Theorems 3.2 and 3.3 in [18] for simplified formulas in the case autonomous state equations 
or with additive noise. 

5.2 Adaptive selection of a time partition 

In order to write a code that automatically determines the time partition (r) h for achieving a prescribed toler- 
ance in the computation of the predictions yt k+1 /t k aud Pt fe+1 /t fc , an adequate adaptive strategy is necessary. 
Since the equations ( 24 )-( 25 1 for the first two conditional moments of y are ordinary differential equations, 



conventional adaptive strategies for numerical integrators of such class of equations are useful. In what follows, 
the adaptive strategy described in [15] is adapted to the LL filter requirements. 
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Once the values for the relative and absolute tolerances rtol y ,rtolp and atol y ,atolp for the first two 
conditional moments, and for the maximum and minimum stepsizes /i max and h m - ln are set, an initial stepsize 
hi needs to be estimated. Specifically, 



where 
with 

and 



*i(u) 



hi = max{h m i a ,mhi{S(y),8(vec(P)),t 1 -t }} 

5(u) = min{10(Wi(u),<5 2 (u)} 
f atol u if do(u) < 10 ■ atol u or di(u) < 10 • atol u 



0.01 



dg(u) 
di(u) 



otherwise 



"max{di(u),d 2 (u)} 
Here, cfo(u) = ||u to / t J, di(u) = \\F(t Q ,u to / to )\\ and tfc(u) 



max{atoZ u , <5i • rtol u }. if max{di{\i), (^(u)} < prs 
0.01 

( r-rrs — , , ) otherwise 



fiF(to,U*o/to) . 9F(t ,u to/to ) 

1 * 1*0) "to /to / 



Ot 



arc 



the norm of the filters and its first two derivatives with respect to t at to, where ||u| 



1 



dim(u) ^ l=1 y st l (vL)> 



with sc 4 (u) = atol u + rtol u 



to/to 



and F is the vector field of the equation for u. 
Starting with the filter estimates Yt k /t k an d ^t k /t k > the basic steps of the adaptive algorithm for determining 
(r) h and computing the predictions yt k+1 /t k an d ~Pt k+1 /t k between two consecutive observations tk and tk+i 
are the following: 

1. Computation of y Tn /t k and vec(P Tn / tk ) at r„ = T n _! + 2h n by the recursive evaluation of the expressions 
( 42 )- ( 43 ) at the two consecutive times r n _i + h n and (r n _i + h n ) + h n . That is, 

yr„/t fc ~ yr„_i/t fc + U T„_i,t fc + u r„_ i+fc n , t fc 

and 

«ec(P t/t J = Lie h » M < T - 1 >u T »- ll t lk + Lie*» M ( T »- 1+h »V B _ 1 +fc B ,iv 

2. Computation of an alternative estimate for the predictions at t„ = r n _i + 2/i„ by means of the expressions 



9r n /t k = yr n -x/t k + L 2 e 



2h„M(r„ 



and 



,(P T „ /t J=L ie 2fc - M (— )u T „_ l!tfc , 



which follow from the straightforward evaluation of ( 42 )-{ 43 1 at r n _i + 2/i„ 
3. Evaluation of the error formulas 



Ex = 



1 



E( 



-) 2 and £ 2 = 



\ d ~{ sci (y) 

where p T „/tk = l ' ec (-f > r„/t fe ) an d sc*(u) = atol u + rio/ u • max{ 



1 f ( ^^ Pr„A M 2 



T„_!/ti 
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4. Estimation of a new stepsize: 



"'new — 

max{/i m i n , mm{S new (Ei),5 new (E2)}} 

where 

!/i„-min{5,max{0.25,0.8-(-)<JTT}} E < 1 
h n • min{l,max{0.1,0.2 • ( — ) w}} £>1 

5. Validation of y Tn /t k and vec(P Tn /t fc ): if max{E\, E%} < 1 or h n = h mirll then accept y r „/t fe an d 
wec(P Tn / tfc ) as approximations to the first two conditional moments of x at r„ = t„_i + 2/i„. Oth- 
erwise, return to step 1 with h n = h new . 

6. Control of the final stepsize: if r„ + 2h n = tk+i, stop. If r„ + 2h n + h new > tfc+i, then redefine 

hnew (j~n ~t~ 2/l n ). 

7. Return to step 1 with n = n + 1 and h n = h new . 

Remarks: It is worth to emphasize that initial stepsize hi is computed just one time for computing the 
value of T\ £ [to, t{\. For other r„ £ [tk, tk+i] with n = n tk + 1 and k > 0, the initial value for the corresponding 
h n is set as h n — h new , where the value h new was estimated when the previous stepsize was accepted. 

Further note that, because the flow property the exponential operator, only two exponential matrices need to 
be evaluated in steps 1 and 2, instead of three. These two exponential matrices can the efficiently computed 
through the well known Pade method for exponential matrices [33 or, alternatively, by means of the Krylov 
subspace method |33j in the case of high dimensional state equation. Even more, low order Pade and Krylov 
methods as suggested in [22] can be used as well for reducing the computation cost, but preserving the order-/3 
of the LL filters. 

5.3 Adaptive LL filter algorithm 

Starting with yt /t = x W*o an d Pt /*o = Qto/to> ^ ne adaptive LL filter is computed by the following recursive 
algorithm of two steps for all k = 0, 1, . . . , N — 1 : 

1. Computation of the predictions y tk+1 /t k £aid Pt k+1 /t k according the adaptive strategy of the previous 
subsection, and evaluating the prediction variance 

~^t k+1 /t k — -ftfc+i/ifc ~ yt k+1 /t k y tk+1 / tk i 

2. Computation of the filters 

yt fc +i/t fc+ i =yt k+ i/t h + K t fc +i( z t fc+1 - c yt k+1 /t k )> 

V tfc+l/tfc + l — V *fc + l/*fc _ K «fc + l CV tfc + l/*fc! 

at tfc+i with filter gain 

K tk+1 = V tfc+l/tfc CT(CV ffc+l/tfc CT + 5VJ" 1 . 
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6 Numerical Simulations 



In this section, the performance of the approximate filters introduced in this paper is illustrated, by means 
of simulations, with four examples of state space models. To do so, the prediction and filtering values are 
computed in four different ways: 1) the exact one, when it is possible; 2) the conventional LL filter; when the 
exact filter is available; 3) the order-1 LL filter with various uniform time partitions; and 4) the adaptive order-1 
LL filter. For each example, the conventional error analysis for the moment estimates and the estimation of 
the weak convergence rate is carried out (see, e.g., [28]. f6]). 

The state space models to be considered are the followings. 

Example 1. State equation with multiplicative noise 

dx = atxdt + o\ftxdw\ (44) 

and observation equation 

z tk =x(t k )+e tk , for k = 0,1,.., N (45) 

with a = —0.1, a = 0.1, t Q = 0.5, £ = 0.0001, x to / tf) = 1 and Qt /t = 0- For this state equation the predictions 
for the first two moments are 

Xt u = x+ u e a( * 2 ~*' )/2 and Q f u = Q t „ e (*+° 2 im 2 -t 2 k ) 

respectively. 

Example 2. State equation with two additive noise 

dx = atxdt + er 1 t p e a * 2/2 G?Wi + a 2 \Ttdw 2 (46) 

and observation equation 

z th =x(t k ) + e th , for A; = 0,1,.., JV (47) 

with a = -0.25, p = 2, a x = 5, a 2 = 0.1, t = 0.01, E = 0.0001, x to/to = 10 and Q to/to = 0. For this state 
equation the predictions for the first two moments are 

T. n - T, u p a (* 2 -*')/ 2 and O, u — (Q* 14 I g 2x a(t 2 -tl) , a \ (f 2p+l f 2p+U at 2 

• L t k+1 /t k — ■ L t k /t k c diiu ^t k+1 /t k — (Wt k /t k + 2p+l 2a' 



respectively. 

Example 3. Van der Pool oscillator with random input [13] 



dxi = x 2 dt (48) 

dx 2 = {~{x\ — l)x 2 — x\ + a)dt + adw (49) 

and observation equation 

z tk =x 1 {t k ) + e tkl for k = 0,1,.., N, (50) 

where a = 0.5 and a 2 = (0.75) 2 are the intensity and the variance of the random input, respectively. In 
addition, t Q = 0, E = 0.001, x tf)/fo = [1, 1], Q t[)/to = 0. 

Example 4. Van der Pool oscillator with random frequency [13] 

dx\ = x 2 dt (51) 

dx 2 — (—(.x 2 — l)x 2 — zoxi)dt + ax\dw (52) 

and observation equation 

zt k = xi(tk) + e tk , for k = 0, 1, ..,N, (53) 
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/i = 1/64 


/i = 1/128 


h = 1/256 


h = 1/512 


P 


h/h 


2.47 ±0.06 x 10" 5 


1.23 ±0.03 x 10" 5 


6.12 ±0.14 x 10 _ti 


3.05 ±0.07 x 10 _ei 


1.01 




8.08 ±0.20 x 10- 6 


4.03 ±0.10 x 10~ y 


2.01 ±0.05 x 10" e 


1.00 ±0.03 x 10~ e 


1.00 


h/H 


4.06 ±0.11 x 10~ e 


2.03 ±0.06 x 10~ y 


1.01 ±0.03 x lO" 6 


5.05 ±0.14 x 10~ 7 


1.00 




2.36 ±0.07 x 10- 6 


1.18 ±0.03 x 10 -6 


5.87 ±0.17 x 10-'' 


2.93 ±0.08 x 10~ 7 


1.00 




1.52 ±0.05 x 10~ H 


7.60 ±0.24 x 10~ 7 


3.78 ±0.12 x 10~ 7 


1.89 ±0.06 x 10~ 7 


1.00 




9.36 ±0.30 x 10" ' 


4.66 ±0.15 x 10~ 7 


2.33 ±0.07 x 10" ' 


1.16 ±0.04 x 10~ 7 


1.00 


tr/tf 


4.49 ±0.22 x 10" 7 


2.23 ±0.10 x 10- 7 


1.11 ±0.05 x 10" 7 


5.55 ±0.27 x 10" 8 


1.00 




1.32 ±0.07 x 10" 7 


6.55 ±0.35 x 10" 8 


3.26 ±0.17 x 10" 8 


1.63 ±0.09 x 10" 8 


1.01 


W*9 


2.42 ±0.19 x 10" 8 


1.19 ±0.09 x 10" 8 


5.94 ±0.46 x 10" 9 


2.96 ± 0.23 x 10" y 


1.01 



Table 1: Confidence interval for the errors between the exact LMV filter x tfc+1 / tfc+1 , U tfc+1( / tfc+1 of ( 44 1- ( 45 I and the 
order-1 LL filter yt k+1 /t k+1 ^t k+1 /t k+1 on ( r )h with different value of h. Order /? of weak convergence estimated from 
the errors. 



where zu = 1 and a 1 = 1 are the frequency mean value and variance, respectively. In addition, to = 0, 
E - 0.001, x to/to = [1, 1] and Q to/to - 0. 

For each example, 2000 realizations of the solution of the state equation were computed by means of 
the Euler [25] or the Local Linearization scheme [55] for the equations with multiplicative or additive noise, 
respectively. For each examples, the realizations where computed over a thin time partition {to + nS : d = 
10~ 4 ,n = 0, ..,9 x 10 4 } for guarantee a precise simulation of the stochastic solutions on the time interval 
[to, to + 9]. A subsample of each realization at the time instants {t}N=9 — {tk = to + k : k = Q,..,N} was 
taken to evaluate the corresponding observation equation. In this way, 2000 time series {zl k }k=o,..,N , with 
« = 1, ..2000, of 10 values each one were finally available for every state space example. 

For each time series of the first two examples, the values of the exact LMV filter, the conventional LL filter 
on {t}pj, the order-1 LL filter on uniform time partition (r) h = {r n = to + nh : n = 0, .., N/h} C {i}./v with 
h = 1/64, 1/128, 1/256, 1/512, and the adaptive order-1 LL filter were computed. 



and 



U t ft " V t It 



and 



IP , -V 1 



For each time series {zl k }k=o...,N, four type of errors were evaluated: the errors 

between each approximate filter and the exact one, and the errors 

between the predictions, for all k = 0, .., N — 1. 
The 2000 errors of each type were arranged into M = 20 batches with K = 100 values each one, which are 
denoted by ejy, / = 1, .., Af; j = l,...,K. Then, the sample mean of the l-th batch and of all batches can be 
computed by 

1 K 1 M 

3=1 1=1 
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i nn 
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h = 1/64 


ft, = 1/128 


h = 1/256 


h = 1/512 


P 


*l/*0 


1.22 ±0.00 x 10" 4 


6.14 ± 0.00 x 10" 5 


3.08 ±0.00 x 10" 5 


1.54 ±0.00 x 10" 5 


1.01 


t 2 /h 


7.61 ± 0.02 x 10~ 5 


3.85 ±0.00 x 10~ b 


1.94 ±0.00 x 10~ 5 


9.71 ±0.02 x 10~ B 


1.00 


h/t 2 


4.05 ±0.04 x 10~ 5 


2.06 ±0.02 x 10~ 5 


1.04 ±0.00 x 10~ 5 


5.22 ±0.05 x 10" y 


1.00 


Ufa 


1.77 ±0.04 x 10~ 5 


9.06 ±0.17 x 10~ B 


4.59 ±0.09 x 10~ B 


2.31 ± 0.04 x 10~ y 


1.00 




6.10 ± 0.14 x 10"° 


3.16 ±0.07 x 10- B 


1.61 ±0.04 x 10- B 


8.09 ±0.18 x 10" 7 


1.00 




1.68 ±0.05 x 10"° 


8.81 ±0.23 x 10~ 7 


4.51 ±0.12 x 10~ 7 


2.28 ±0.06 x 10~ 7 


1.00 


h/h 


3.79 ±0.11 x 10- 7 


1.99 ±0.06 x 10" 7 


1.02 ±0.03 x 10~ 7 


5.17 ±0.15 x 10" 8 


1.00 


h/tr 


8.01 ±0.34 x 10" 8 


4.14 ±0.18 x 10" 8 


2.10 ±0.09 x 10" 8 


1.06 ±0.04 x 10" 8 


1.00 


W*8 


1.58 ±0.07 x 10" 8 


7.83 ±0.35 x 10" y 


3.90 ±0.17 x 10" y 


1.95 ±0.09 x 10" y 


1.00 



Table 2: Confidence interval for the errors between the exact LMV predictions x tfc+1 / t 



, U tfe+l/tfc of (|44f-(|45} and their 
approximations yt k+1 /t k ^t k+1 /t k obtained by the order-1 LL filter on (r)^ with different value of h. Order f3 of weak 
convergence estimated from the errors. 



respectively. The confidence interval for each type of error is computed as 

[e - A, e + A], 

where 



1 M 



M e M - 1 ^ 1 



A = tl-a/2,Af-l\/ VT: = — — r > \et — e 



1=1 



and ii- Q /2,M-i denotes the 1 — a/2 percentile of the Student's t distribution with M — 1 degrees for the 
significance level < a < 1. The 90% confidence interval (i.e., the values A for a — 0.1) was chosen. 



6.1 Results for Example 1 



Tables T]|3 show the estimated errors for the state space model (44 1- ([45]). Specifically, Table [I] shows the 



confidence interval for the errors between the exact LMV filter x t u ,XJ t u and the order-1 LL 
filter y tfc+l/tfc+1 , V tfc+l/tfc+1 on the time partition (r)£, with h = 1/64,1/128,1/256,1/512. Table [| shows the 
confidence interval for the errors between the exact LMV predictions x tfc+1 / tfc , V tk+1 / tk and their approximations 
yt k+1 /t k 7 ^t k+1 /t k obtained by the order-1 LL filter on {t)\. Table [3] shows the confidence interval for the errors 
between the moments of the exact LMV filter and their respective approximations obtained by the conventional 
LL filter and the adaptive LL filter. The number of accepted and fail steps of the adaptive LL filter at each 
tk € {t}pf is given in Figure 1. The absolute and relative tolerances for the first and second moments were set 
as rtol y = rtolp = 5 x 10~ 9 and atol y = 5 x 10~ 9 , atolp — 5 x 10~ 12 . Note as the accuracy of the LL filter on 
uniform partitions (r)^ improve as h decrease, and the large difference among the accuracy of the conventional 
and the adaptive LL filter. 
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k 


n 
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A 

V / 


tk+i/tk 


■\tA 
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2.79 ±0.00 x 10~ 3 
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5.62 ± 0.13 x 10 J 
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2.09 ± 0.00 x 10 ° 


2 
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4.04 ± 0.07 x 10 


1.16 ± 0.01 x 10 


3 


5.74 ± 0.06 x 10 a 


6.75 ± 0.12 x 10 8 
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r- c\r\ i f\ -i t\ i r\ — 7 

5.29 ± 0.10 x 10 ' 


4 


4.54 ± 0.05 x 10 


5.73 ± 0.17 x 10 5 


1.70 ± 0.04 x 10 


1.92 ± 0.04 x 10 


5 


3.17 ± 0.04 x 10 J 


5.72 ± 0.15 x 10 8 


7.21 ± 0.20 x 10 


5.62 ± 0.15 x 10 


6 


2.01 ± 0.03 x 10 


6.07 ± 0.14 x 10 


2.44 ± 0.07 x 10 


1.34 ± 0.04 x 10 


7 


1.24 ± 0.02 x 10 ^ 


5.02 ± 0.15 x 10 8 


7.37 ± 0.27 x 10 


a or 1 A 1 A -i /A — Q 

2.85 ± 0.12 x 10 9 


8 


7.32 ± 0.14 x 10 


2.82 ± 0.07 x 10 8 


1.81 ± 0.09 x 10 


r- aa i r\r»o -\ r\ — 111 

5.39 ± 0.23 x 10 u 


k 


o 

V / 

J tfc+i/tfe+i 


A 

•^tfc+l/ifc+l 


V y 


XT A 

*fc+i/*fc+i 





3.94 ±0.08 x 10~ 3 


3.50 ±0.07 x 10~ 7 


7.22 ±0.17 x 10~ 3 


6.38 ±0.14 x 10~ 7 


1 


6.25 ±0.13 x 10~ 4 


1.43 ±0.03 x 10~ 7 


9.61 ±0.25 x 10~ 4 


2.16 ±0.05 x 10~ 7 


2 


3.58 ±0.07 x 10" 4 


1.01 ±0.02 x 10~ 7 


4.12 ±0.12 x 10~ 4 


1.12 ±0.03 x 10" 7 


3 


3.09 ±0.09 x 10~ 4 


9.44 ±0.29 x 10~ 8 


2.33 ±0.07 x 10~ 4 


6.69 ±0.19 x 10~ 8 


4 


3.50 ±0.09 x 10~ 4 


1.04 ±0.03 x 10~ 7 


1.64 ± 0.06 x 10" 4 


4.45 ±0.14 x 10" 8 


5 


4.49 ±0.12 x 10~ 4 


1.22 ±0.03 x 10~ 7 


1.16 ±0.04 x 10~ 4 


2.86 ±0.09 x 10~ 8 


6 


5.93 ±0.10 x 10" 4 


1.12 ±0.03 x 10" 7 


7.96 ±0.28 x 10" 5 


1.45 ±0.07 x 10~ 8 


7 


6.61 ±0.12 x 10" 4 


6.93 ±0.18 x 10" 8 


3.92 ±0.17 x 10" 5 


4.50 ±0.23 x 10" 9 


8 


5.91 ±0.09 x 10" 4 


2.89 ±0.11 x 10~ 8 


1.48 ±0.06 x 10" 5 


8.17 ±0.63 x 10" 1U 



Table 3: Confidence interval for the errors between the exact LMV filter and predictions of (44l-(45l with thier 
corresponding approximations obtained by the conventional LL filter and the adaptative LL filter, which are denoted 
with superscripts and A, respectively. 



For each approximate conditional moment, the estimated order /? of weak convergence were obtained as the 
slope of the straight line fitted to the set of points {log 2 (hj), log 2 (e(/ij))} j=1 4 taken from their corresponding 

errors Tables and § The values j3 are shown in these tables as well. The estimates /3 ~ 1 corroborate 
the theoretical value for f3 given in Theorem 10 The estimate (3 = 2.00 corresponding to y^/ to m Table 
[2] agrees with the expected estimated of j3 for the equation (44) on [£o,*i]< In this particular situation, the 
exact prediction x tl / to given by (|3j) reduces to an ordinary differential equation and the LL prediction formula 
24l reduces to the classical order-2 LL integrator for such class of equations (see, e.g., [20 ). In the others 



subintervals [i/t.i/c+i] with k 7^ 0, the prediction yf fc+1 / t( . depends nonlinearly of y through the initial value 



k+i/tk+i 



6.2 Results for Example 2 



Tables 4](6 show the estimated errors for the state space model ([46])-( 47 ) . In particular, Table [4] shows the 
confidence interval for the errors between the exact LMV filter x t / t , Ut k+1 /t k+1 and the order-1 LL 
filter y* fc+1 /t fc+1 , V tfc+1 / tfc+1 on the time partition {t)\, with h — 1/64,1/128,1/256,1/512. Table [5] shows the 
confidence interval for the errors between the exact LMV predictions x tfc+1 / tfc , V tk+1 / tk and their approximations 



ytfc+i/tic i Vt k+1 /t k obtained by the order-1 LL filter on (r)^. Table [6] shows the confidence interval for the errors 
between the moments of the exact LMV filter and their respective approximations obtained by the conventional 
LL filter and the adaptive LL filter. The number of accepted and fail steps of the adaptive LL filter at each 
tk G {Ojv is given in Figure 1. The absolute and relative tolerances for the first and second moments for this 
filter were set as rtol y — rtolp — 5 x 10~ 8 and atol y — 5 x 10 -8 , atolp = 5 x 10~ n . Note as the accuracy of 
the LL filter on uniform partitions {r) h improve as h decrease, and the large difference among the accuracy of 
the conventional and the adaptive LL filter. 
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h/H 
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1.73 ±0.08 x 10-'' 
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tr/tf 
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1.85 ±0.07 x 10- 7 
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0.97 




2.35 ±0.13 x 10" 7 


1.21 ± 0.07 x 10" 7 


6.11 ± 0.34 x 10" 8 


3.08 ±0.17 x 10" 8 


0.98 


W*9 


1.67 ±0.09 x 10- 7 


8.31 ±0.04 x 10" 8 
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1.00 



Table 4: Confidence interval for the errors between the exact LMV filter x tfc+1 / tfc+1 , U tfc+1 / tfc+1 of (46i-(47l and the 
order-1 LL filter yt k+1 /t k+1 ^t k+1 /t k+1 on ( r )h with different value of h. Order /? of weak convergence estimated from 
the errors. 



For each approximate conditional moment, the estimated order /3 of weak convergence were obtained as the 
slope of the straight line fitted to the set of points {\og 2 (hj), log 2 (e(/ij))} j=1 4 taken from their corresponding 

errors Tables [i] and ^ The values j3 are included in these tables too. The estimates /3 ~ 1 corroborate the 

The estimate (3 ss 2.00 corresponding to y 1 } in Table 5 agrees 
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theoretical value for /? given in Theorem 
with the expected estimated of /3 for the > 

the exact prediction x tfe+1 / tfc given by p|) reduces to an ordinary differential equation and the LL prediction 



equation (46) on [tktk+i], for all k. Similarly to the previous example, 



formula (24) reduces as well to the classical order-2 LL integrator for all k. Contrary to the first example, in 
this one, the prediction y^ k+1 / tk with k ^ does not depend of y through the initial value yf fc+1 / tfc+1 and so 

the estimate j3 « 2.00 is preserved. 



6.3 Results for Examples 3 and 4 



Since explicit formulas of the LMV filter for the state space models (48)-(50) and ( 51 )-( 53 ) are not available, 
the error analysis of the previous examples should be adjusted. In this situation, by taking into account the 
results of the previous examples, the moments estimated by the adaptive LL filter with small tolerance can be 
used as a precise estimation for the moments of the exact LMV filter. By doing so, the confidence interval for 
the errors can similarly be computed as before for estimate the order f3 of weak convergence of the order-1 LL 
filter. Table [7] shows the estimated order /3 of weak convergence obtained, as explained above, as the slope of 
the straight line fitted to the set of points {\og 2 (hj), log 2 (e(/ij))} j=1 4 , where e(hj) denotes the error between 
order-1 LL filter on (r)^, with h = 1/64, 1/128, 1/256, 1/512, and the adaptive LL filter with small tolerance. 
The tolerances for adaptive filter were set as rtol y = rtolp = 5 x 10~ 8 and atol y = 5 x 10~ 8 , atolp = 5 x 10~ n 
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h = 1/64 


h = 1/128 


h = 1/256 


h = 1/512 


P 


*l/*0 


2.43 ±0.00 x 10 -3 


1.28 ±0.00 x 10" 3 


6.56 ±0.00 x 10" 4 


3.32 ±0.00 x 10" 4 


0.88 


t 2 /h 


7.22 ±0.00 x 10~ 2 


3.54 ±0.00 x 10~ 2 


1.75 ±0.00 x 10~ 2 


8.73 ±0.00 x 10- 3 


1.01 


h/t 2 


1.69 ±0.00 x 10- 1 


8.29 ±0.00 x 10~ 2 


4.11 ±0.00 x 10~ 2 


2.04 ±0.00 x 10~ 2 


1.01 


Ufa 


1.16 ±0.00 x 10" 1 


5.73 ±0.00 x 10~ 2 


2.84 ±0.00 x 10~ 2 


1.42 ±0.00 x 10~ 2 


1.01 




3.38 ±0.00 x 10~ 2 


1.68 ±0.00 x 10~ 2 


8.36 ±0.00 x 10- a 


4.17 ±0.00 x lO" 3 


1.00 




4.81 ±0.00 x 10~ 3 


2.41 ± 0.00 x 10" 3 


1.21 ±0.00 x 10" 3 


6.05 ±0.00 x 10~ 4 


0.99 


h/h 


3.77 ±0.00 x 10" 4 


1.91 ±0.00 x 10" 4 


9.62 ±0.00 x 10~ 5 


4.83 ±0.00 x 10- 5 


0.97 


h/tr 


3.27 ±0.00 x 10- 5 


1.65 ±0.00 x 10" 5 


8.28 ±0.00 x 10" 6 


4.15 ±0.00 x 10" 6 


0.98 


W*8 


1.70 ± 0.00 x 10" 5 


8.44 ± 0.00 x lO" 6 


4.21 ± 0.00 x 10" B 


2.10 ±0.00 x 10" 6 


1.00 



Table 5: Confidence interval for the errors between the exact LMV predictions x tfc+1 / t 



, U t(c+l/tfc of (|46f-(|47f and their 
approximations Yt k+1 /t k > ^t k+1 /t k obtained by the order-1 LL filter on (r)^ with different value of h. Order f3 of weak 
convergence estimated from the errors. 



for the model ( 48 1-( 50 ) , and as rtol y — rtolp = 1 x 10 7 and atol y = 1 x 10 7 , atolp = 1 x 10 10 for the model 
(51 1- (53 1. For each model, the number of accepted and fail steps of the adaptive LL filter at each tk G {t}N is 
given in Figure 1. Notice that, for both examples, the estimates (3 « 1 corroborate the theoretical value for /3 
stated in Theorem [TOl 



7 Conclusions 

Approximate Linear Minimum Variance filters for continuous-discrete state space models were introduced and 
their order of convergence stated. As particular instance, the order-/3 Local Linearization filters were studied in 
detail. For them, practical algorithms were also provided and their performance in simulation illustrated with 
various examples. Simulations shown that: 1) on thin time partitions between observations, the order-1 LL 
filter provides satisfactory approximations to the exact LMV filter; 2) the convergence of the approximate LL 
filter to the exact LMV filter when the maximum stepsize of time partition between observations decrease; 3) 
the adaptive LL filter provides an automatic and computationally efficient approximation to the LMV filtering 
problem; and 4) the effectiveness of the approximate LL filter for the identification of nonlinear stochastic 
system from a reduced number of partial and noisy observations. 
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k 


y + , . , /+, 

v L j^ + I ] L 


z k+l/ T k 


V° , 

c fc + l / z k 


c fc + l/ l k 





7.69 ±0.00 x 10~ 2 


2.17 ± 0.00 x 10~ 6 


2.63 ±0.00 


3.72 ± 0.00 x 10" 4 


1 


2.09 ±0.01 x 10 _1 


2.14 ±0.04 x 10" 7 


8.01 ±0.03 


1.85 ± 0.00 x 10" 3 


2 


2.81 ± 0.06 x 10 _1 


8.41 ±0.38 x 10~ 8 


4.93 ±0.13 x 10 2 


3.24 ± 0.02 x 10" 3 


3 


4.02 ±0.07 x 10" 1 


1.26 ±0.07 x 10~ 7 


3.22 ±0.17 x 10 2 


2.33 ± 0.02 x 10" 3 


4 


3.82 ±0.08 x 10 _1 


1.55 ± 0.08 x 10~ 7 


6.18 ±0.10 x 10 1 


7.45 ± 0.06 x 10" 4 


5 


2.27 ± 0.06 x 10" 1 


1.06 ±0.06 x 10" 7 


6.23 ±0.25 x 10" 1 


1.15 ± 0.01 x 10" 4 


6 


9.36 ± 0.21 x 10" 2 


4.68 ±0.22 x 10~ 8 


8.34 ±0.25 x 10" 2 


9.74 ± 0.11 x 10" 6 


7 


2.89 ±0.07 x 10" 2 


1.23 ±0.08 x 10~ 8 


8.19 ±0.28 x 10" 3 


7.43 ± 0.08 x 10-'' 


8 


8.63 ±0.28 x 10~ 3 


1.10 ± 0.06 x 10" 9 


2.73 ±0.02 x 10~ 3 


2.86 ±0.00 x 10~ 7 


k 


v° , 

•> ifc+i/tfc-fi 


v A , 

J tk + l/tk + 1 


V° , 

tk + l/tk + 1 


t-k + lft-k + 1 





1.75 ±0.03 x 10~ 3 


3.29 ±0.06 x 10~ 9 


3.08 ±0.08 x 10- 2 


5.73 ±0.14 x 10" 8 


1 


9.47 ±0.24 x 10" 7 


3.37 ±0.07 x 10" 1U 


1.62 ±0.07 x 10" 5 


6.84 ±0.26 x 10" 9 


2 


2.17 ±0.06 x 10" 6 


2.16 ±0.05 x 10" 10 


5.58 ±0.23 x 10" 5 


5.75 ±0.25 x 10" 9 


3 


2.73 ±0.04 x 10" 6 


3.15 ±0.05 x 10~ 1U 


5.70 ±0.17 x 10" b 


7.10 ±0.26 x 10" 9 


4 


3.02 ±0.05 x 10" 6 


6.58 ±0.14 x 10^ 10 


2.79 ±0.11 x KT 5 


8.01 ±0.42 x 10" 9 


5 


1.35 ±0.03 x 10" b 


1.94 ± 0.06 x 10~ 9 


7.37 ±0.28 x 10" b 


9.32 ±0.44 x 10" 9 


6 


1.16 ±0.14 x 10" 3 


7.07 ±0.20 x 10" 9 


1.89 ±0.09 x KT 4 


1.04 ±0.06 x 10" 8 


7 


1.01 ± 0.03 x 10" 3 


1.26 ±0.05 x 10" 8 


3.18 ±0.16 x KT 5 


5.56 ±0.36 x 10" 9 


8 


6.91 ±0.18 x 10" 5 


8.05 ±0.21 x 10" 9 


2.01 ±0.12 x KT 5 


2.83 ±0.16 x 10" 9 



Table 6: Confidence interval for the errors between the exact LMV filter and predictions of 1 46 1-1 47 1 with thier 
corresponding approximations obtained by the conventional LL filter and the adaptative LL filter, which are denoted 
with superscripts and A, respectivelly. 



k\Add 




v *fc+i/tfc 


y*A; + l/tA; + l 


+ 1 /tk + 1 


k\Mul 


ytk+i/tk 


v tfc+i/tfc 


ytk+lf tk+1 







1.10 


1.04 


1.11 


1.04 





1.08 


1.01 


1.03 


1.01 


1 


1.04 


1.05 


1.04 


1.05 


1 


1.03 


1.03 


1.04 


1.03 


2 


1.03 


1.03 


1.03 


1.03 


2 


1.03 


1.04 


1.05 


1.04 


3 


1.02 


1.02 


1.03 


1.02 


3 


1.03 


1.03 


1.04 


1.02 


4 


1.01 


1.01 


1.01 


0.97 


4 


1.02 


1.02 


1.02 


1.02 


5 


1.01 


1.03 


1.01 


1.01 


5 


1.02 


0.94 


1.02 


0.83 


6 


1.02 


1.01 


1.01 


0.98 


6 


1.01 


0.98 


1.01 


0.97 


7 


1.02 


1.04 


1.02 


1.06 


7 


0.97 


1.00 


1.02 


0.99 


8 


1.03 


1.02 


1.02 


1.02 


8 


1.02 


1.01 


1.03 


0.99 



Table 7: Estimate order of convergence P for the moments of the order-1 LL filter applied to the state space models 
( |48[ )-(|50| and (|51[)-(|53[) corresponding to the Van der Pool oscillator with additive (Add) and multiplicative (Mul) noise, 
respectivelly. 
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Figure 1: Number of accepted and fail steps of the adaptive LL filter at each tk € {t}N for each one of the 
four examples. 
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